📖 第16课：Lattice Gas自动机

LatticeGas HPP 流体力学 Navier-Stokes

🌊 用CA模拟流体——Lattice Gas的诞生

1986年，Frisch、Hasslacher和Pomeau发表了划时代的论文——他们证明了可以用极其简单的CA来模拟流体力学！这就是Lattice Gas自动机（LGA），它从微观粒子碰撞出发，在宏观尺度上再现了Navier-Stokes方程。

核心思想：

流体 = 大量粒子的集体行为

每个粒子在网格上运动，碰撞时遵循守恒律

宏观密度和速度 = 微观粒子的统计平均

当粒子数足够多时 → 恢复Navier-Stokes方程！

HPP模型——最早的Lattice Gas

HPP模型（Hardy, Pomeau, de Pazzis, 1973）：

正方形网格，Von Neumann邻域
4个速度方向：东、南、西、北
每个格点最多4个粒子（每个方向至多1个）
碰撞规则：对头碰撞的两个粒子旋转90°

碰撞规则（唯一非平凡情况）：

东西方向各有一个粒子 → 变为南北方向各有一个粒子

南北方向各有一个粒子 → 变为东西方向各有一个粒子

守恒律：粒子数守恒 + 动量守恒

HPP模型实现

// ============================================================================
// lattice_gas_hpp.v - HPP Lattice Gas自动机
// 正方形网格，4速度方向，弹性碰撞
// ============================================================================
module lattice_gas_hpp #(
    parameter WIDTH  = 64,
    parameter HEIGHT = 64
)(
    input  wire             clk,
    input  wire             rst_n,
    input  wire             enable,
    input  wire             init,
    // 4位/格点：位0=东, 位1=南, 位2=西, 位3=北
    output wire [WIDTH*HEIGHT*4-1:0] state,
    output wire [31:0]      step_count,
    output wire [31:0]      total_particles
);

    // 每格点4位状态：[北,西,南,东]
    reg [3:0] grid [0:WIDTH*HEIGHT-1];
    reg [3:0] grid_nxt [0:WIDTH*HEIGHT-1];
    reg [31:0] steps;

    integer idx, x, y;

    always @(*) begin
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
            x = idx % WIDTH;
            y = idx / WIDTH;

            // ---- 步骤1：碰撞（在同一格点内） ----
            reg [3:0] after_collision;
            case (grid[idx])
                4'b0101: after_collision = 4'b1010;  // 东西→南北
                4'b1010: after_collision = 4'b0101;  // 南北→东西
                default: after_collision = grid[idx];  // 其他不变
            endcase

            // ---- 步骤2：传播（移动到邻居格点） ----
            // 东向粒子(x) → 移到(x+1,y)
            // 南向粒子(y) → 移到(x,y+1)
            // 西向粒子(x) → 移到(x-1,y)
            // 北向粒子(y) → 移到(x,y-1)
            // 即：从邻居接收粒子

            // 来自西邻居的东向粒子
            wire from_west  = (x > 0) ? grid[y*WIDTH+x-1][0] : grid[y*WIDTH+WIDTH-1][0];
            // 来自北邻居的南向粒子
            wire from_north = (y > 0) ? grid[(y-1)*WIDTH+x][1] : grid[(HEIGHT-1)*WIDTH+x][1];
            // 来自东邻居的西向粒子
            wire from_east  = (x < WIDTH-1) ? grid[y*WIDTH+x+1][2] : grid[y*WIDTH][2];
            // 来自南邻居的北向粒子
            wire from_south = (y < HEIGHT-1) ? grid[(y+1)*WIDTH+x][3] : grid[x][3];

            grid_nxt[idx] = {from_south, from_east, from_north, from_west};
        end
    end

    // 状态更新
    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) begin
            steps <= 32'd0;
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                grid[idx] <= 4'd0;
        end else if (init) begin
            steps <= 32'd0;
            // 初始配置：随机粒子
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                grid[idx] <= 4'd0;  // 空白
        end else if (enable) begin
            for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
                grid[idx] <= grid_nxt[idx];
            steps <= steps + 32'd1;
        end
    end

    // 粒子总数
    reg [31:0] total;
    always @(*) begin
        total = 0;
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1)
            total = total + grid[idx][0] + grid[idx][1] +
                    grid[idx][2] + grid[idx][3];
    end

    // 输出
    genvar gi;
    generate
        for (gi = 0; gi < WIDTH*HEIGHT; gi = gi + 1)
            assign state[gi*4 +: 4] = grid[gi];
    endgenerate
    assign step_count = steps;
    assign total_particles = total;

endmodule

HPP的缺陷

HPP模型有一个根本缺陷——缺乏旋转不变性。正方形网格导致流体在0°/90°方向的行为与45°方向不同，这在宏观上表现为非物理的各向异性。

HPP的宏观方程：

密度：ρ = n₀ + n₁ + n₂ + n₃（4方向粒子数之和）

动量：ρv = (n₀-n₂)x̂ + (n₁-n₃)ŷ

问题：应力张量不是各向同性的，导致宏观方程偏离Navier-Stokes

💡 解决方案：FHP模型（下一课）使用三角形网格（六角形邻域），解决了各向异性问题。这是Lattice Gas理论的重大突破。

🏋️ 练习

练习16.1：在HPP模型中创建一个"风洞"——左边界持续注入东向粒子，右边界排出。观察粒子流的稳态模式。

练习16.2：在HPP网格中放置一个"障碍物"（禁止粒子进入的格点），观察流动如何绕过障碍。是否形成类似卡门涡街的图案？

练习16.3：验证HPP的守恒律——每步后粒子总数和总动量是否保持不变？

练习16.4：修改碰撞规则，使其包含"三体碰撞"。这对宏观行为有什么影响？

练习16.5（挑战）：从HPP的微观规则推导宏观的动量方程。证明应力张量是各向异性的。

🏆 成就解锁

🏅 流体仿真先驱

你已经进入物理仿真阶段！HPP模型虽然简单，但包含了Lattice Gas的核心思想——从微观碰撞到宏观流体的涌现。

✅ 实现HPP Lattice Gas引擎
✅ 理解碰撞-传播的两步过程
✅ 掌握守恒律在CA中的编码
✅ 了解HPP的各向异性缺陷

📐 深入分析：Lattice Gas的守恒律验证

守恒律的数学表述：

粒子数守恒：Σ_{all cells} (n₀+n₁+n₂+n₃) = const

动量守恒：Σ_{all cells} [(n₀-n₂)x̂ + (n₁-n₃)ŷ] = const

这两个守恒律是Lattice Gas宏观行为的基础——它们保证了宏观的质量和动量方程

守恒律的硬件验证：

每步后计算全网格的粒子总数和总动量

与上一步比较，任何差异都表示实现错误

这是最有效的debug工具——如果守恒律被违反，必然有bug

补充实现

// 守恒律验证器 - 实时验证Lattice Gas的守恒律
module conservation_checker #(
    parameter WIDTH = 64,
    parameter HEIGHT = 64
)(
    input  wire             clk,
    input  wire             rst_n,
    input  wire [3:0]       grid [0:WIDTH*HEIGHT-1],
    output wire [31:0]      total_particles,
    output wire signed [31:0] momentum_x,
    output wire signed [31:0] momentum_y,
    output wire             conservation_ok  // 1=守恒律满足
);
    reg [31:0] particles;
    reg signed [31:0] px, py;
    reg [31:0] prev_particles;
    reg signed [31:0] prev_px, prev_py;
    integer idx;
    
    always @(posedge clk or negedge rst_n) begin
        if (!rst_n) begin
            prev_particles <= 32'd0;
            prev_px <= 32'd0; prev_py <= 32'd0;
        end else begin
            prev_particles <= particles;
            prev_px <= px; prev_py <= py;
        end
    end
    
    always @(*) begin
        particles = 32'd0; px = 32'd0; py = 32'd0;
        for (idx = 0; idx < WIDTH*HEIGHT; idx = idx + 1) begin
            particles = particles + grid[idx][0] + grid[idx][1] +
                        grid[idx][2] + grid[idx][3];
            px = px + grid[idx][0] - grid[idx][2];  // 东-西
            py = py + grid[idx][1] - grid[idx][3];  // 南-北
        end
    end
    
    assign total_particles = particles;
    assign momentum_x = px;
    assign momentum_y = py;
    assign conservation_ok = (particles == prev_particles);
endmodule

性能与优化分析

HPP vs FHP 的定量对比：

正方形网格(HPP)：应力张量有4个独立分量，不满足旋转不变性

三角形网格(FHP)：应力张量满足旋转不变性，等价于2个独立分量

宏观差异：HPP的粘度张量是各向异性的，FHP的粘度是标量

这意味着FHP可以正确模拟圆管中的Poiseuille流，而HPP不行

📖 扩展阅读

Wolfram, S. A New Kind of Science, Wolfram Media, 2002.
Toffoli, T. & Margolus, N. Cellular Automata Machines, MIT Press, 1987.
Ilachinski, A. Cellular Automata: A Discrete Universe, World Scientific, 2001.

🔬 补充专题：实验方法论

在进行CA实验时，科学的方法论至关重要。以下是一些通用的实验指导原则：

实验设计三要素：

控制变量：每次只改变一个参数，保持其他不变
可重复性：记录所有参数，确保实验可以精确复现
统计显著性：多次运行取平均，避免偶然结果误导

参数扫描方法

系统性地扫描参数空间是理解CA行为的关键技术：

参数	范围	步长	测量指标
规则号	0-255	1	种群密度/周期
初始密度	0.1-0.9	0.1	收敛时间
网格大小	16-256	×2	有限尺寸效应
边界条件	环形/固定/镜像	离散	边界效应

💡 实验记录模板：每次实验记录以下信息——日期、规则号、初始条件、网格大小、边界条件、运行步数、关键观察、数据文件路径。这是科学研究的良好习惯，也使得结果可以追溯和验证。

数据可视化技巧

CA实验产生的数据通常是高维的（空间+时间+状态）。有效的可视化对于理解至关重要：

时空图：横轴=空间位置，纵轴=时间，颜色=状态——最常用的CA可视化方式
相图：x轴=参数A，y轴=参数B，颜色=测量指标——用于参数空间扫描
频谱图：对种群时间序列做FFT，识别周期性
直方图：统计状态/密度的分布

自相关分析：

时间自相关：R(τ) = ⟨n(t)·n(t+τ)⟩ / ⟨n²⟩

空间自相关：C(r) = ⟨n(x)·n(x+r)⟩ / ⟨n²⟩

如果R(τ)以周期T振荡 → 系统有周期T的行为

如果C(r)幂律衰减 → 系统处于临界状态

◀ 第15课 📜 目录第17课 ▶

元胞自动机课程 · 从Conway到Langton到Lattice Gas